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Abstract: A complete error calibration model with 12 independent parameters is established 
by analyzing the three-axis magnetometer error mechanism. The said model conforms to an 
ellipsoid restriction, the parameters of the ellipsoid equation are estimated, and the ellipsoid 
coefficient matrix is derived. However, the calibration matrix cannot be determined 
completely, as there are fewer ellipsoid parameters than calibration model parameters. 
Mathematically, the calibration matrix derived from the ellipsoid coefficient matrix by a 
different matrix decomposition method is not unique, and there exists an unknown rotation 
matrix R between them. This paper puts forward a constant intersection angle method 
(angles between the geomagnetic field and gravitational field are fixed) to estimate R. The 
Tikhonov method is adopted to solve the problem that rounding errors or other errors may 
seriously affect the calculation results of R when the condition number of the matrix is very 
large. The geomagnetic field vector and heading error are further corrected by R. The 
constant intersection angle method is convenient and practical, as it is free from any 
additional calibration procedure or coordinate transformation. In addition, the simulation 
experiment indicates that the heading error declines from +1° calibrated by classical 
ellipsoid fitting to +0.2° calibrated by a constant intersection angle method, and the 
signal-to-noise ratio is 50 dB. The actual experiment exhibits that the heading error is further 
corrected from +0.8° calibrated by the classical ellipsoid fitting to +0.3° calibrated by a 
constant intersection angle method. 
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1. Introduction 

Accurate measurement of the geomagnetic field has been widely applied in geophysical research, 
space magnetic measurement, military defense, mineral resources exploration, and drilling practice [1]. 
Common magnetometers include a helium optical-pumping magnetometer, a proton magnetometer, a 
three-axis flux gate magnetometer, and an anisotropic magnetoresistive magnetometer. The former two 
magnetometers could only be used to measure geomagnetic field magnitude, whereas the latter two 
magnetometers, belonging to the vector magnetometer type, are the most common magnetic sensors in 
an attitude and heading reference system (AHRS), which is characterized by low cost and high 
reliability. However, the magnetometer is full of inevitable errors such as zero deviation, sensitivity 
errors, non-orthogonal errors, misalignment errors, measurement noise, hard iron errors, and soft iron 
errors, and consequently there are large errors between the measurement results and actual geomagnetic 
field vector. As a result, calibration and error compensation should be conducted before the application. 
The calibration of magnetometers can be divided into the geomagnetic domain calibration and the 
heading domain calibration. The heading domain calibration only calibrates the heading error of the 
magnetic compass composed of the magnetometer, but each axis output of the magnetometer has not 
been calibrated. The geomagnetic domain calibration directly calibrates the output of each axis of the 
magnetometer. The magnetometer output after calibration can not only calculate the accurate heading 
angle (a = tm''^(H" / H")) but also can be used in geomagnetic navigation for providing rich 
component information. 

The "two-step" algorithm directly calibrates the magnetometer error without an external heading 
reference. The principle of the two-step algorithm is that when the magnitude of the true geomagnetic 
field is invariable at the calibration spot, the ideal locus of the three-axis magnetometer output lies on a 
spherical surface, whereas the actual locus of the magnetometer measurement outcome lies on an 
ellipsoid under the disturbance of error. The first step of the two-step algorithm transforms the 
non-linear observation equation based on the square of the magnetic field magnitude into a linear 
equation of the new unknown parameter by a mathematical transformation and acquires the unknown 
parameter by a batch least squares solution. In the second step, the biases and the scale factors can be 
derived from an inverse analytical solution [2], but the non-orthogonal and misalignment errors cannot 
be calibrated because of the error model, which includes only zero deviation and sensitivity error. The 
iterative batch least squares method is used in [3] to obtain the ellipsoid parameter unbiased estimation, 
from which the initial conditions can be acquired by the aforementioned two-step algorithm. The 
algorithm possesses consistent convergence, and the posterior covariance could be used as a metric for 
the compensation effect. The "extension of the two-step" algorithm improves the two-step calibration 
algorithm by compensating for the non-orthogonal error [4]. However, as the assumed sensor x-axis 
coincides with the body x-axis (it is also possible that none of the three axes coincides in reality), the 
misalignment error cannot be calibrated. 
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It is proposed in literature [5] that there are only nine independent parameters in the fitting ellipsoid 
model. Failing to clarify the 12 parameters in the calibration model mathematically could also be 
interpreted to mean that there is more than one method to decompose the ellipsoid coefficient matrix into 
two orthogonal matrixes (error parameter matriK). According to matrix theory, an orthogonal rotation 
matrix (indicating the rotation transform in three-dimensional space) exists between matrixes, resulting 
from different decomposition methods. In other words, magnetometer data can be compensated for 
sensor errors and the presence of magnetic distortions by mapping an ellipsoid of data to a sphere; 
however, the rotation of the sphere is unknown [5,6]. The essence of this difference discloses the failure 
to clarify the three-axis misalignment error. It has been shown in [7] that the magnetometer 
measurement noise is zero-mean Gaussian white noise, but the noise of the square of the magnetometer 
output is not zero mean, and ordinary least square (OLS) estimation is an inconsistent biased estimation. 
Adaptive least square (ALS) is utilized to obtain a consistent unbiased estimator of ellipsoid parameters 
because of the continuous adjustment procedure to the OLS cost function. This adjustment requires a 
known noise variance. Consistent estimation to an unknown noise variance, however, can be realized [8]. 
ALS improves the accuracy of the ellipsoid fitting, but the calibration parameters derived by eigen 
decomposition are not always true because of the different matrix decomposition methods mentioned 
above. In addition, the coefficient matrix of the ellipsoid fitting could also be used to derive the error 
parameter matrix from singular value decomposition [9]. However, the error model of this method 
approximates a symmetric matrix, which contains only nine undetermined coefficients; thus, the 
application of this method has limitations. In [10] a geomagnetic field vector and some auxiliary vector 
dot products are choosen as the constants to exam the effect of the rotation matrix in rectifying and 
compensating the error. This method does not require additional calibration, and the compensation 
accuracy depends on the choice of auxiliary vector. When the auxiliary vector is nearly parallel with the 
geomagnetic vector, the worst compensation result is obtained, whereas the best result is obtained when 
the former vector is perpendicular to the latter vector. Because the auxiliary vector is a constant vector of 
the earth coordinate system, it should be converted into components of the body coordinate system in a 
mathematical operation. Maximum likelihood estimation is used in the error parameter fitting, and 
misalignment error requires additional independent calibration for identification [11]. 

The transition matrix defined in [12] between the sensor coordinate system and body system contains 
six unknown parameters. This matrix completely describes the non-orthogonal error and misalignment 
error and adopts the "three-step" algorithm in magnetometer pre-calibration before flight, without 
reference to any direction information. The misalignment rotation matrix can be derived from the 
condition that the vertical component of the geomagnetic field remains unchanged while rotating the 
magnetometer respectively around the x-axis, y-axis, and z-axis in the geoid. It is proposed in [13] that 
the rotation matrix of the misalignment error could be obtained when the sensor rotates around an axis, 
and the projection of the magnetic field output is invariable, which lies on the plane that is perpendicular 
to the said axis. This method is more valuable in engineering as it is unnecessary to keep leveling. This 
paper is based on the constant intersection angle between the geomagnetic field vector and gravitation 
field vector [14], without additional independent calibration and coordinate system conversion, and 
seeks to derive a rotation matrix more conveniently, furthering the compensation of magnetic field 
measurement error and heading error. 
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2. Error Modeling 

Because of the influence of the manufacturing processes and environment, a three-axis magnetometer 
has manufacturing, installation and environmental errors. The manufacturing error contains a 
zero-deviation Hq and sensitivity error installation error includes non-orthogonal error /Cnor ^iid 
misalignment error ; environmental error includes the soft iron error and hard iron error //p. On 
account of all the errors listed above, the relationship between the true geomagnetic field the 
three-axis magnetometer output is demonstrated in Equation (1), in which indicates the Gaussian 
white noise is zero-mean and its standard deviation is a: 

The different types of errors are specifically shown in Subsections 2. 1-2.6. 
2.1. Sensitivity Error 

Sensitivity represents the proportional relationship between the magnetometer input and output. 
Sensitivity error results from the inconsistency between the amplification of the measuring electronic 
circuit and its nominal value, which can be expressed as: 

K^=dmg{K^K^yKJ (2) 

If sensitivity is the only error that affects the sensor, the Equation (2) refers to the measurement 
H^ = K H 



e 



2.2. Zero Deviation 

There is often a small voltage bias in the sensor output signal, even when there is no magnetic 
intensity. In such conditions, the sensor produces non-zero output, which can be expressed as: 

Ho=[H,, Hoy Hj (3) 
If zero deviation is the only error that affects the sensor, the Equation (3) refers to the measurement 

Hm — He + Hq. 

2.3. Non-Orthogonal Error 

Non-orthogonal error is caused by the limitation of machining, which leads the x-, y- and z-axes to not 
being completely orthogonal to each other. This is shown in Figure 1, assuming that the z-axis of the 
sensor completely coincides with the z-axis of the orthogonal coordinate system 7^: 

cos a 0 sin a 

sinpcosy cosPcosy siny (4) 
0 0 1 



If non-orthogonal error is the only error that affects the sensor. Equation (4) refers to the 
measurement i/„ = K„ ■ . 
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Figure 1. Sensor coordinate system non-orthogonal error. 

iiZm 




Xm 

2.4. Misalignment Error 

During the installation, the sensor x-axis may not totally coincide with the body vertical axis, as 
illustrated in Figure 2. The resultant transition matrix (rotation matrix) 7^ exists between the virtual 
orthogonal coordinate system 7^ and the body orthogonal coordinate system 7^, which is equivalent to 
a slight angle generated from rotations of the sensor around the x-axis, y-axis and z-axis, respectively: 





ixy 


^XZ 


iyx 




ty2 


izx 


tzy 


^ZZ 



(5) 



If misalignment is the only error that affects the sensor, the Equation (2) refers to the measurement 



Figure 2. Body orthogonal coordinate system and virtual orthogonal coordinate system. 

Zm ♦ Zb 





Xm 



2.5. Soft Iron Errors 

Ferromagnetic materials produce an inductive magnetic field under the influence of a geomagnetic 
field and the induced electric current around the sensor. It also varies with the attitude of the 
magnetometer and the position in the magnetic field, which can be expressed as: 





'0-xx 


0-xy 


^xz' 




ClyX 


^yy 


Uyz 




P'ZX 


0-zy 


^zz. 



(6) 
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The Qij tenns represent the effective soft iron coefficients and are the constants of proportionality 
between the magnetic field applied to soft iron and the resulting induced magnetic field [3]. For example, 
Uxy represents the effective coefficient relating the field generated in the x-direction in response to an 
applied field in the y-direction, where the effective soft iron coefficients appears insignificant when 
i ^ j. Mathematically, the soft iron error for principal diagonal elements is equivalent to the sensitivity 
error; and non-diagonal elements correspond to the non-orthogonal error and misalignment error. If soft 
iron is the only error that affects the sensor. Equation (2) takes to the measurement //m = Cs • H^. 

2.6. Hard Iron Errors 

This type of error results from permanent magnets and magnetic hysteresis, that is, remanence of 
magnetized iron materials. Mathematically, this error is equivalent to zero-deviation: 

//p = [V V ^v^V (7) 
If hard iron is the only error that affects the sensor. Equation (2) refers to the measurement 

Ignoring the measurement noise in Equation (1), the geomagnetic field vector can be derived as: 

He = G-iH^- Ho) - Hp (8) 

dll 9l2 913' 



in which G - M - 



921 922 923 
931 932 933 



3. Classical Ellipsoid Assumption Algorithm 

When an ideal magnetometer rotates in the calibration spot, the geomagnetic field magnitude remains 
constant, and the output locus of the three-axis magnetometer is a sphere. In contrast, a non-ideal 
magnetometer output locus would be an ellipsoid because of the influence of error. The square of H^ is 
derived from Equation (8), and sorted into the quadratic form of the surface: 



,T G'G ^ ^ nyC + ^G ^ UyGH, + 2^GH, + U'^H^ 

"fef " \\H\f " \\H\f 



(9) 



Equation (9) is rewritten as a general form of the quadratic surface equation: 

H^^H^) = ^Hl+bH^H^^+cHi^+dH^H^+eH^^H^ + JHl^+pH^+qH^^+rH^+s = 0 (lO) 

where ^ = (a, b, c, d, e, f, p, q, r, sY. 

Generally, non-diagonal elements of the soft iron matrix , the non- orthogonal error and 
misalignment error of the soft iron error, are minuscle, therefore, the error matrix M is strictly diagonally 
dominant [5]. Therefore, G^G is also strictly diagonally dominant, belonging to a positive-definite real 
symmetric matrix, and Equation (9) is a quadratic ellipsoid equation, then for convenience. Equation (9) 
can be written as: 

X^AX-2X^AX + XlAXo = l (11) 



Sensors 2014, 14 



8491 



is the coefficient matrix after ellipsoid fitting, and Xq = —A ^ 



are 



a b/2 dll 
where, A= b/2 c ell 
dll ell f 

the spherical center coordinates. The objective of magnetometer calibration is to conduct ellipsoid fitting 
based on the set of measurement data, deriving the coefficient matrix A and spherical center 
coordinates Xq. 

The principle of ellipsoid fitting is finding an ideal ellipsoid where the sum of squares of the algebraic 
distances calculated from the measured data point to the ellipsoid is minimized [15], which is expressed 
as Equation (12) or (13): 



N 



mm 



in^||F(^H^i)lP 



1=1 



where: D = 



mxl^^myl 



min t' D'^DE 



(12) 



(13) 



'mzl "mxl 



mzl 



ij2 



In order to assure that the fitting result is an ellipsoid, Nikos contends that parameter o should 
guarantee that matrix A in Equation (10) is positive definite matrix or negative matrix [16]. This 
constraint is analyzed in literature [16], which states that the equality constraint 4ab — b^ = 1 can be 
imposed for the convenience of calculation, which would not lose its generality and can be expressed in 
the matrix form: 



where C = 



^3x7 



and Ci = 



0 



fC^ = 1 



(14) 



By introducing a Lagrange multiplier X, we can obtain the simultaneous equation from Equation (13) 
and Equation (14) as follows: 

G(^) = ^^D'^D^ + X(l - a'^CQ (15) 

In order to obtain the extremum under constraints, the first-order derivative could be derived from 
Equation (15), and Equation (16) could be obtained when the first-order derivative of Equation (15) is 
set to zero: 

D^D^ = XC^ (16) 
Therefore, the extremum under an ellipsoid constraint could be obtained: 

I f = 1 

It has been proved in [17] that only one of the generalized eigenvalues in Equation (16) is positive, 
and its eigenvector is the best ellipsoid fitting parameter ^. The ellipsoid coefficient matrix A and 
spherical center coordinate Xq is derived from further calculation. It could be concluded from the 
comparison between Equation (9) and Equation (10) that: 



(17) 
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G'^G = llHelP 

Hp = G(Zo-Ho) (1^) 

||//el| is not necessarily the actual magnitude of the geomagnetic field because its magnitude only 
affects the components on the x-axis and y-axis, whereas the heading hinges on the proportion of the 
geomagnetic field projection on the x-coordinate and y-coordinate [9]. Introducing =G(Xq -H^) 

to Equation (8), we obtain H^=G(H^-H^)-G(X^-H^) = G(H^-X^) . Therefore, it IS not 

necessary to use Hq. The 3x3 matrix G combining scale factors, non-orthogonal errors, misalignments, 
and soft iron disturbances is not strict symmetrical, including nine independent parameters, and there are 
three combined biases in the error model. The total number of error parameters should be 12. 

Therefore, the key to deriving calibration matrix G is to decompose the ellipsoid coefficient matrix 
into two orthogonal matrixes. However, the calibration matrix cannot be clarified, as the orthogonal 
matrix resulting from a different decomposition is not exactly the same. Because G is also a positive 
definite matrix, the decomposition results Gl and G2 only differ in an orthogonal rotation matrix R, 
which is equivalent to a fixed rotating angle in the body coordinate system between the geomagnetic 
field vector after compensation and the true geomagnetic field vector [5]. 

This paper provides a constant intersection angle method to solve the rotation matrix R. On the basis 
of this rotation matrix, the matrix G resulting from orthogonal decomposition could be transformed into 
the calibration matrix that is needed. The calibrated geomagnetic field could be obtained by solving 
Equation (8) and Equation (18) simultaneously: 

He = Q-CH^-Xo) (19) 

where Q = R - G. 

4. Constant Intersection Angle Method 

Similar to the classical ellipsoid fitting method, the constant intersection angle method also requires 
calibration in the same location in order to satisfy the condition that the magnitude and direction of the 
geomagnetic field and gravitation field are invariant. The intersection angle between the geomagnetic 
field vector and gravitation field vector at the same measurement spot is constant [14], and the dot 
product of magnetic field vector He and gravity acceleration vector Ag, which ||//g||-y4g -cosG , is 

invariable and can be calculated from the model of geomagnetic field and gravitation field. In the said 
equation, 6 is the intersection angle between the geomagnetic field vector and gravitation field vector, 
and its value can be interpreted to have approximately the same effect as WH^W in Equation (18). 
Therefore, different values of 6 would only affect the magnitude of the geomagnetic field after 
calibration rather than the heading solution. In addition, the acceleration coordinate system is also the body 
coordinate system as the accelerometer is fixed on the body and keeps moving with the body. Thus the 
intersection angle between the geomagnetic field vector and gravitation field vector Ag at the same 

measurement spot is constant, even if there is no transformation into a geodetic coordinate system. 
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4.1. Rotation Matrix Solution 



To begin, calibration matrix G is obtained by a singular value decomposition of the positive matrix 
||//ell^ ■ ^ generated from Equation (18). Then, can be derived as below, by introducing G andZp 
into Equation (19) and preliminarily calibrating the geomagnetic field: 

H, = G-CH^-Xo) (20) 



where the gravity acceleration vector is Ag = 



1 



, the magnetic field vector is He = R • H^, and R as 



the undetermined rotation matrix can be represented asR = 



^11 ^12^13 
^21^22^23 
^31^32^33 



Additionally, the dot product 



of the magnetic field vector and the gravity acceleration vector could is described as: 





[A 1 


T 


^11 ^12^13 








A 




^21^22^23 




Hf.y 




Agz. 




-^31^32^33- 







= \\He\\ ■ lUsll ■ cos6=const 



(21) 



Recasting Equation (21) as a linear combination of a changing input vector and an estimation 
parameter vector. Equation (22) is obtained: 



const 



1 [/ 

^gx^cy 

^gx"cz 

A H 

^gy"cx 

A H 

A U 

^gy"cz 

A H 

"gz''cx 

^gz^cy 
A H 



Rii 

^12 
^13 
^21 
^22 
^23 
^31 
^32 
^33 



(22) 



The matrix containing As for the set of A'^ input vectors Ag and H^, is expressed as Equation (23): 



rconsti 



const 



A M 

"■gxl^^cxl 

A hi 

•"gxl^'cyl 

A M 
A M 

A M 

^gyl^'cyl 

A H 
•"gyl^'czl 

A M 
•"gzl^'cxl 

A M 

"^gzl"cyl 

Agzl^czl 



A H 

"^gxn"cxn 

A H 
^gxn'^cyn 

A H 

^gxn^^czn 

A H 

gyn^ cxn 

A H 

^gyn''cyn 

A H 

"^gyn"czn 

A H 

^gzn"cxn 

A H 

•"gzn^'cyn 

A H 
^gzn^^czn 



Rll 
^12 
^13 
^21 
^22 
^23 
^31 
^32 
^33 



(23) 



(p) ^ • (p^ -Y is derived by batch a least squares solution. In the 



Y =(p-x. 

The least squares solution x = (cp^ 
calculation process according to the method proposed in [6] we find the problem that the condition 
number of matrix (p is larger, which decrease the estimation accuracy of the rotation matrix. Therefore, 
the Tikhonov method is adopted to obtain a better estimator. 
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4.2. Noise Suppression 

A high rotation matrix accuracy plays an important role in further correcting the heading error, which 
is calibrated by the classical ellipsoid fitting method. The measurement noise, however, are inevitable 
and have influence on the estimator of rotation matrix R. The noise suppression methodology can reduce 
the sensitivity of the solution to measurement noise. 

When the condition number of matrix (p in Equation (23) is very large, the matrix is approximately 
singular, or when observation data is inaccurate, the solution will be submerged by errors. The Tikhonov 
regularization method could solve the problem by adding other information about the solution, so as to 
identify a meaningful, stable solution [18]. The principle of adding information is that the 2-nonn of the 
solution should be minimized. The regular solution defined by the regularization method is shown 
below [18]: 

mm||(p;c-r||i + X'||Lx||i (24) 

Where L is generally a unit matrix, is a regularization parameter that controls the sensitivity of 9, 
and Y is the disturbance caused by the solution. The selection of 1 mainly relies on experience to 
guarantee matrix is of full column rank. A least squares solution is derived as the solution of 

Equation (24) is equivalent to the solution of Equation (25): 



Equation (25) can be derived from the QR (decompose a matrix into a orthogonal and an upper triangular 
matrix) orthogonal least square algorithm, which is equivalent to solving the following Unear equations: 

((P^(P + ^Vl> = (P^F (26) 

USV^ is then obtained by singular value decomposition of the coefficient matrix (<p^<p + X^L^L). 
The coefficient matrix is a square matrix; matrix U is an orthogonal matrix and matrix 5 is a diagonal 
matrix containing the eigenvalues, V= U. Then x = (USV^ ^ cp^j = V'^i' 'f/ 'cp^j and rotation matrix 

R will obtain by recasting matrix x . Though the noise suppression can reduce the influence of 
measurement noise, other methods are needed to ensure that R is close to a rotation matrix. For example, 
a quite large number of samples can improve the estimation accuracy. Another method is set a constraint 
R^R = I on matrix R , and try to minimize to obtain rotation matrix R. where R^is recast by 

matrix x ,R h the constrained optimization rotation matrix, its determinant is one. R is more close to the 
true rotation matrix i?o than R^- Introducing it into Equation (19) or the equation Hg = R • Hf., heading 
error will be calibrated completely. 

5. Simulation Study 

During magnetometer calibration, to ensure the quality of the measurement data set, a proper uniform 
distribution of the data over all directions must be accomplished [19]. Firstly, 46 data points are 
collected for the simulation, which are evenly distributed on the sphere. According to the principles of 
the equivalent intersection angle, regarding the sphere center as the coordinate origin, and the radius of 
the sphere as one, the x-axis, y-axis and z-axis coordinates are collected in the D46raw matrix, which is 
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composed of 46 sample points. The collection method is shown in Figure 3, where the intersection 
angles AD1OD2 and /.D^OD-^ between any two adjacent data points are equal. Secondly, 46 data points 
are collected for the acceleration data, of which the magnitude is 9.8. As demonstrated in Figure 4, 
assuming that the acceleration vector and unit sphere intersects each other inA^ andAj, ADiOAi is 
equivalent to /LDjOAj. 

Figure 3. Collecting data according to the equivalent intersection angle principle. 




Figure 4. Constant intersection angle between the geomagnetic field vector and gravitation field vector. 

Zffi 




Ym 



Then assuming that error matrix M = 



, zero deviation Hq = 



0.03 
0.02 
0.06 



, and hard 



iron error Hp — 



the output of the three-axis magnetometer influenced by the errors is 



1.06 0.02 0.03 
0.01 0.98 0.04 
0.02 0.05 1.1 

0.0651 
0.02 
0.05 

raw ■ which is distributed on the disturbing ellipsoid. Then, the Gaussian white 
noise of A^(0, 0.002 ) is added into £'46enipse' where 0.002 is the standard deviation of the noise. 

After determining the ellipsoid coefficient matrix, the preliminary calibration matrix G is derived by 
singular value decomposition. Regarding the heading as a reference, which is calculated from the 
original data set D46raw that is free from the influence of the error matrix, the heading error calibrated 
by matrix G is shown in Figure 5a. The heading error is ±1° where the blue dashed Une refers to the 



D46 



ellipse 



D46 
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heading error curve after a preliminary calibration of matrix G, which is derived from the classical 
ellipsoid fitting. Because the ellipsoid model after fitting only contains nine independent parameters, it 
can't determine the 12 error parameters in the error model completely. In other words, the classical 
ellipsoid fitting method can transform the ellipsoid into a sphere, where transformed sphere and the 
original sphere differ by an undetermined rotation matrix R. Consequently, this influences the accuracy 
of the heading calibration. The proposed constant intersection angle method can derive R. 

The constant intersection angle method is not necessary for the additional calibration procedure 
mentioned in reference [11,12] and the coordinate transformation mentioned in reference [10]. Rotation 
matrix R can be derived from three-axis magnetometer and accelerometer data. The calibration matrix G 
derived from ellipsoid fitting coefficients is further corrected by rotation R, allowing us to improve the 
estimation of the calibration matrix Q. the heading error calibrated by matrix Q is lower than +0.2°, 
which smaller than the heading error of +1° when calibrated by classical ellipsoid fitting, shown as the 
red solid line in Figure 5a. 

Figure 5. (a) Heading error when the standard deviation of the noise is 0.002; (b) Heading 
error when the standard deviation of the noise is 0.003; and (c) Heading error when the 
standard deviation of the noise is 0.005. 
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Figure 5. Cont. 
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When the standard deviation of the Gaussian white noise is 0.003, and the magnitude of the 
simulation data is one, the heading error shown is shown in Figure 5b. The heading error calibrated by 
the constant intersection angle method decreases to +0.4° from the value of +1° when calibrated by the 
classical ellipsoid fitting. When the standard deviation of the noise is 0.005, the heading error is also 
shown in Figure 5c. The heading error decreases from +1.2° to +0.6° when calibrated by the constant 
intersection angle method. The error increases with the noise. 

6. Experiments 

6.1. Measurement Equipment and Measurement Spot 

In order to compare the calibration results between the classical ellipsoid fitting method and the 
constant intersection angle method proposed earlier, the nonmagnetic turntable experiment adopts a 
self-developed magnetic compass, which is mounted on the nonmagnetic turntable. The magnetic 
compass adopts a PIC24HJ128GP504 microprocessor (Microchip Technology, Chandler AZ, USA), 
equipped with a fluxgate sensor as the magnetic sensing element, and a Honeywell QA-T160 
accelerometer (Honeywell Conglomerate Company, Morristown, NJ, USA) is used. The output heading 
angle of the optical-electrical encoder is regarded as the true heading reference, where the error is less 
than 3'. As shown in Figure 6, the magnetic compass is fixed on the nonmagnetic turntable. Where the 
body coordinate frame {b) takes the forward direction in the magnetic compass as the direction of Xb-axis, 
taking the downward direction as the Zb-axis, and setting the jb-axis in conformity with the left-hand rule, 
the navigation coordinate frame (n) takes a direction in the north-east-downward frame. It can be 
inferred from the body coordinate frame and navigation coordinate frame that the clockwise heading 
angle is positive, and the heading angle is 0° when the JCb-axis points northward. 

The experiment was performed in a Beijing suburb, where the magnetic environment in the vicinity 
of the magnetometer is better than that in the laboratory in the city. The principle of the constant 
intersection angle is adopted to ensure that the data points are distributed evenly on the sphere. The data 
collection is composed of 46 sets of data spread on the x-axis, };-axis and z-axis, collected as follows. 
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To begin, where the magnetic compass is fixed on the nonmagnetic turntable, and the roll angle is 0°, 
and pitch angle is 60° and -60°, respectively, the nonmagnetic turntable is rotated clockwise such that 
the magnetic compass JCb-axis can intersect with the north direction at angles of 6°, 66°, 126°, 186°, 246° 
and 306°; Where the roll angle is 0°, and the pitch angle is respectively 60° and -60°, the nonmagnetic 
turntable is rotated clockwise such that the magnetometer JCb- axis can intersect with the north direction at 
angles of 6°, 42°, 78°, 1 14°, 150°, 186°, 222°, 258°, 294° and 330°. Where both the roll angle and the 
pitch angle are 0°, the nonmagnetic turntable is rotated clockwise such that the magnetometer Xb-axis can 
intersect with the north direction at angles of 6°, 36°, 66°, 96°, 126°, 156°, 186°, 216°, 246°, 276°, 306°, 
and 336°. Where the roll angle is 0°, and the pitch angle is 90° and -90°, the .x:b-axis is vertically upward 
and downward, respectively. The accelerometer data should be collected simultaneously at the same 
46 positions, because the method is attitude-dependent. 

Figure 6. Nonmagnetic turntable. 




Then, magnetometer preliminary calibration can be conducted through the ellipsoid assumption 
method. The heading angle \\i could be is derived Equation (27). It should be noted that magnetic field 
data in the body coordinate frame {H^, Hy, H^) should be transformed to the navigation coordinate 
frame (H^, Hy, H^). The transformation matrix is shown in Equation (28): 

'a {H^> 0, > 0) 

n + a (H^ <0) 

2n + a {H^ <0, HJ}> O) 



¥ = 



{H^ >0,H^ = 0) 



71 

2 

y(//]^<0, = 0) 



(27) 



where a = tan'^ (Hi^/H^). 
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(28) 



6 and y respectively refer to the pitch angle and roll angle, which are determined by the three-axis 
accelerations gx, Qy, and g^, when either angle is not zero and where g = ^g^ + gy + gz- 

Q = -sm-^{3-/g) 



= -sin-(^V,) 



(29) 



(30) 



6.2. Accelerometer Calibration 



As a type of inertial sensor, accelerometers are reliable, independent, and immune to environmental 
impacts. The primary errors of accelerometers are bias, sensitivity error, and non-orthogonal error. Such 
errors are similar to the magnetometer errors, as observed in Equation (1-3). The following equation 
describes the accelerometer observation: 

gm = K-Kor-gQ+go = C-gQ+gQ (31) 

where, is the bias of accelerometer, g^ is the actual acceleration of gravity, and 

Sy 0 aSr 

is derived from Equation (1) and Equation (3). Generally, the non-orthogonal 



C 



0 



0 



angle of a three-axis accelerometer provided by the manufacturer is less than 0.5°. For calculation 
purposes, the trigonometric functions of Equation (3) can be simplified. The actual acceleration of 
gravity is derived from Equation (31) as follows: 

-C-^iSm-go) (32) 



The output trajectory of the three-axis accelerometer accords with the ellipsoid hypothesis as well. 
There are nine unknown parameters consisting of three sensitivity errors, three non-orthogonal error and 
three bias. Thus, the error parameters of the accelerometer can all be derived from the coefficients of the 
ellipsoid equation. The experiment is performed on the turntable shown in Figure 6 by mounting the 
magnetic compass on the aligned turntable, aligning the body frame of the magnetic compass with the 
turntable by using the location pins, and then rotating the turntable to different angular positions. The 
results in Figure 7 show that the absolute error declines from 0.1582 m/s to 4.6323 x e m/s , and the 
mean square root error drops from 0.0545 m/s to 1.0358 x e m/s ; thus, the maximum error of the 
accelerometer can be reduced by a factor of approximately 1,000. 
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Figure 7. The accelerometer error before calibration and after calibration. 
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6.3. Experimental Result 

The calibration matrix G of Equation (18) is derived from the ellipsoid coefficient matrix by 
eigen-value decomposition and singular value decomposition in the classical ellipsoid fitting method. 
Then G is introduced into Equation (20) to calibrate the measurement of the geomagnetic field and 
obtain H^. calibrated by the classical ellipsoid fitting method. However, calibration matrix G is not the 
actual calibration matrix, which should be multiplied by the rotation matrix R mentioned earlier. R is 
obtained using the constant intersection angle method to further improve the accuracy of magnetic 
sensor. This method can be used to other type of three-axis inertia sensor including misalignment error. 
Then, the complete calibration geomagnetic field = R ■ H^. is obtained. The heading error is shown in 
Figure 8 a. The red solid line refers to the heading error calibrated by the constant intersection angle 
method, which is less than +0.3°, and the blue dashed line refers to the heading error calibrated by the 
classical ellipsoid fitting method, which is within the range of +0.8°. 

Figures 8b,c are the results of two other experiments conducted in the laboratory in the city, where the 
red solid line refers to the heading error calibrated by the constant intersection angle method, and the 
blue dashed line refers to the heading error calibrated by the classical ellipsoid fitting method. The 
maximum heading errors calibrated by the constant intersection angle method are +0.35° in Figure 8b 
and +0.4° in Figure 8c. The maximum heading errors calibrated by the classical ellipsoid fitting method 
are approximately +0.9° in Figure 8b and +1.1° in Figure 8c. This shows that the constant intersection angle 
method has similar calibration performance in the laboratory in the city. The heading errors after calibration 
are slightly larger because magnetic perturbations are larger in the city. 

As mentioned in Section 4.2, the accuracy of rotation matrix R affects the heading error directly. The 
measurement noises, however, are inevitable and have influence on the estimator of rotation matrix R. 
The noise suppression methodology can reduce the sensitivity of the solution to measurement noise, and 
the constrained optimization (set a constraint R^R = I on matrix R) can further improve the estimation 
accuracy of R. In order to verify the effect of this methodology, we conduct a comparison experiment 
shown in Figure 9. 

The green line refers to the heading error that the rotation matrix R is estimated directly by the 
measurement data; the red line is the heading error that the estimation of R is adopted the noise suppression 
methodology. It can be seen from Figure 9 that the heading error decrease from +0.48° to +0.32°. 
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Figure 8. (a) Experiment results comparing heading error in Beijing suburbs; 
(b) experimental results comparing heading error in laboratory in the city; and (c) more 
experimental results comparing heading error in the laboratory in the city. 
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Figure 9. Comparison experiment about the noise suppression methodology. 
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7. Conclusions 

This study thoroughly analyses the three-axis magnetometer measurement errors and establishes a 
complete error model including 12 independent parameters, which is more universal and conforms to the 
ellipsoid restriction. However, the calibration matrix derived from the ellipsoid coefficient matrix by a 
different matrix decomposition method is not unique, and there exists an unknown rotation matrix R 
between them [5]. A constant intersection angle method (angles between geomagnetic vector and 
gravitational vector are fixed in the same location) is proposed to estimate R. The geomagnetic field 
vector and heading error are further corrected by R. The simulation experiment indicates that the heading 
error declines from +1° when calibrated by the classical ellipsoid fitting method to +0.2° when calibrated 
by the constant intersection angle method, and the signal-to-noise-ratio is 50 dB. The actual experiment 
demonstrates that the heading error is further corrected by the constant intersection angle method 
decreases from the value of ±0.8° when calibrated by the classical ellipsoid fitting method to a value of 
+0.3°. The method proposed is convenient and practical, as it is free from any additional calibration 
procedure mentioned in references [11,12] and the coordinate transformation mentioned in reference [10]. 
In addition, the noise suppression methodology can reduce the sensitivity of the solution to measurement 
noise and obtain higher accurate estimator of the rotation matrix than reference [6]. 

Acknowledgments 

This paper is supported by "Projects from National Natural Science Foundation of China" 
(No.61273082) and "Fundamental Research Funds for the Central Universities" (No.FRF-TP-09-017B). 
Furthermore, the authors would also like to thank all who have helped with this study. 

Author Contributions 

Yan Xia Liu writes most of the paper. Xi Sheng Li designs the structure of the paper and help to 
revise the manuscript. Xiao Juan Zhang provides the experiment data. Yi Bo Feng corrects part of the 
English expression. All authors approve the final version of the manuscript. 



Sensors 2014, 14 



8503 



Conflicts of Interest 

The authors declare no conflict of interest. 
References 

1 . Zhang, Z. ; Cheng, D. ; Lian, M. ; Zhou, Z. ; Wang, J. Design of the optically pumped magnetometer' s 
signal detection and control loop. Electron. Meas. lustrum. 2011, 25, 366-371. 

2. Gebre-Egziabher, D.; Elkaim, G.H.; Powell, J.D.; Parkinson, B.W. A non-linear, two-step 
estimation algorithm for calibrating solid-state strapdown magnetometers. In Proceedings of the 8th 
International Conference on Navigation Systems, St. Petersburg, Russia, 27-31 May 2001; 
pp. 200-299. 

3. Gebre-Egziabher, D.; Elkaim, G.H. Calibration of strapdown magnetometers in magnetic field 
domain. /. Aerosp. Eng. 2006, 19, 87-102. 

4. Foster, C.C.; Elkaim, G.H. Extension of a two-step calibration methodology to include 
nonorthogonal sensor axes. IEEE Trans. Aerosp. Electron. Syst. 2008, 44, 1070-1078. 

5. Li, Z.; Li, X. Research on magnetic deviation compensation of three-axis electronic compass based 
on ellipsoid hypothesis. Chin. J. Sci. Instrum. 2011, 32, 2210-2215. 

6. Kok, M.; Hoi, J.D.; Schon, T.B.; Gustafsson, F.; Luinge, H. Calibration of a magnetometer in 
combination with inertial sensors. In Proceedings of the 15th International Conference on 
Information Fusion (FUSION), Singapore, Singapore, 9-12 July 2012; pp. 787-793. 

7. Renaudin, V.; Afzl, M.H.; Lachapelle, G. Complete triaxis magnetometer calibration in the 
magnetic domain. J. Sens. 2010, 2010, 1-10. 

8. Kukush, A.; Markovsky, I.; van Huffel, S. Consistent estimation in an implicit quadratic 
measurement error model. Comput. Stat. Data Anal. 2004, 47, \13-\A1. 

9. Fang, J.C.; Sun, H.W.; Cao, J.; Zhang, X.; Tao. Y. A novel calibration method of magnetic compass 
based on ellipsoid fitting. IEEE Trans. Instrum. Meas. 2011, 60, 2053-2061. 

10. Li, X.; Li, Z. Dot product invariance method for the calibration of three-axis magnetometer in 
attitude and heading reference system. Chin. J. Sci. Instrum. 2012, 33, 1813-1818. 

11. Vasconcelos, J. F.; Elkaim, G.; Silvestre, C; Oliveira, P.; Cardeira, B. Geometric approach to 
strapdown magnetometer calibration in sensor frame. IEEE Trans. Aerosp. Electron. Syst. 2011, 47, 
1293-1306. 

12. Asl, H.G.; Pourtakdoust, S.H.; Samani, M. A new non-linear algorithm for complete pre-flight 
calibration of magnetometers in the geomagnetic field domain. Proc. Inst. Mech. Eng. G J. Aerosp. Eng. 
2009, 223, 729-739. 

13. Bonnet, S.; Bassompierre, C; Godin, C; Lesecq, S.; Barraud, A. Calibration methods for inertial 
and magnetic sensors. Sens. Actuators A Phys. 2009, 156, 302-311. 

14. Sodhi, R.; Prunty, J.; Hsu, G.; Oh, B. Automatic Calibration of A Three-Axis Magnetic Compass. 
U.S. Patent 7,451,549 Bl, 18 November 2008. 

15. Feng, W.; Liu, S.B.; Liu, S.W.; Yang, S. A calibration method of three-axis magnetic sensor based 
on ellipsoid fitting. /. Inf. Comput. Sci. 2013, 10, 1551-1558. 



Sensors 2014, 14 



8504 



16. Nikos, G.; Michael, G.S. Head detection and tracking by 2-D and 3-D ellipsoid fitting. 
In Proceedings of the International Conference on Computer Graphics, Geneva, Switzerland, 
19-24 June 2000; pp. 221-226. 

17. Halir, R.; Flusser, J. Numerically stable direct least squares fitting of ellipses. 

Proc. Int. Conf. Cent. Eur. Comput. Gr. 1998, 1, 125—132. 

18. Wang, J.-G.; Wei. T.; Zhou, Y.-B. Tikhonov regularization method for a backward problem for the 
time-fractional diffusion equation. Appl. Math. Model. 2013, 37, 8518-8532. 

19. Merayo, J.M.G.; Brauer, P.; Primdahl, F.; Petersen, J.R.; Nielsen, O.V. Scalar calibration of vector 
magnetometers. Meas. Sci. Technol. 2000, 11, 120-132. 



© 2014 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article 
distributed under the terms and conditions of the Creative Commons Attribution license 
(http://creativecommons.Org/licenses/by/3.0/). 



